{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2.1.2 Building blocks for programming preconditioners"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import netgen.gui\n",
    "from netgen.geom2d import unit_square\n",
    "from ngsolve import *\n",
    "%gui tk\n",
    "from ngsolve.la import EigenValues_Preconditioner\n",
    "mesh = Mesh(unit_square.GenerateMesh(maxh=0.1))\n",
    "Draw (mesh)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "fes = H1(mesh, order=3, dirichlet=\"left|bottom\")\n",
    "u, v = fes.TnT()\n",
    "a = BilinearForm(fes)\n",
    "a += SymbolicBFI(grad(u)*grad(v) + u*v)\n",
    "a.Assemble()\n",
    "f = LinearForm(fes)\n",
    "f += SymbolicLFI(1*v)\n",
    "f.Assemble()\n",
    "gfu = GridFunction(fes)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Create a Jacobi-preconditioner\n",
    "\n",
    "Let  $A=$ `a.mat` be the assembled matrix, which can be decomposed based on `FreeDofs` ($F$) the remainder ($D$), as in $\\S$[1.3](../unit-1.3-dirichlet/dirichlet.ipynb),\n",
    "\n",
    "$$\n",
    "A = \\left( \\begin{array}{cc} A_{FF} & A_{FD} \\\\ A_{DF} & A_{DD} \\end{array} \\right). \n",
    "$$\n",
    "\n",
    "Then the matrix form of the **point Jacobi preconditioner** is\n",
    "\n",
    "$$\n",
    "J = \\left( \\begin{array}{cc} \\text{diag}(A_{FF})^{-1} & 0  \\\\ 0 & I  \\end{array} \\right),\n",
    "$$\n",
    "\n",
    "which can be obtained in NGSolve using `CreateSmoother`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "preJpoint = a.mat.CreateSmoother(fes.FreeDofs())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "NGSolve also gives us a facility to quickly check an estimate of the condition number of the preconditioned matrix by applying the Lanczos algorithm on the preconditioned system.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       " 0.0146975\n",
       " 0.0567568\n",
       " 0.0954973\n",
       " 0.112599\n",
       " 0.149013\n",
       " 0.198016\n",
       " 0.266185\n",
       " 0.354506\n",
       " 0.452415\n",
       " 0.535142\n",
       " 0.63906\n",
       " 0.751009\n",
       " 0.836039\n",
       " 0.932691\n",
       " 1.04077\n",
       " 1.16473\n",
       " 1.24579\n",
       " 1.36685\n",
       " 1.47433\n",
       " 1.58784\n",
       " 1.70789\n",
       " 1.83511\n",
       " 1.92562\n",
       " 2.02092\n",
       " 2.12651\n",
       " 2.20727\n",
       " 2.27824\n",
       " 2.38882\n",
       " 2.44463\n",
       " 2.47644\n",
       " 2.51031\n",
       " 2.54019\n",
       "  2.8097"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lams = EigenValues_Preconditioner(mat=a.mat, pre=preJpoint)\n",
    "lams"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "An estimate of the condition number\n",
    "\n",
    "$$\n",
    "\\kappa = \\frac{\\lambda_{\\text{max}} }{ \\lambda_{\\text{min}}}\n",
    "$$\n",
    "\n",
    "is therefore given as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "191.168278526202"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "max(lams)/min(lams)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One might wonder if we have gained anything by this point Jacobi preconditioning. What if we did not precondition at all?\n",
    "\n",
    "Not preconditioning is the same as preconditioning by the identity operator on $F$-dofs. One way to realize this identity operator in NGSolve is through the projection into the space of free dofs (i.e., the space spanned by the shape functions corresponding to free dofs). NGSolve provides\n",
    "\n",
    "```\n",
    "Projector(mask, range)   # mask: bit array, range: bool\n",
    "```\n",
    "\n",
    "which projects into the space spanned by the shape functions of the degrees of freedom marked as range in mask."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "preI = Projector(mask=fes.FreeDofs(), range=True)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "*Note*  that another way to obtain the identity matrix in NGSolve is    \n",
    "```\n",
    "IdentityMatrix(fes.ndof, complex=False).\n",
    "```"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "1585.7752877045514"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lams = EigenValues_Preconditioner(mat=a.mat, pre=preI)\n",
    "max(lams)/min(lams)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Clearly the point Jacobi preconditioner has reduced the condition number. \n",
    "\n",
    "We can use preconditioners within iterative solvers provided by NGSolve's `solvers` module (which has `MinRes`, `QMR` etc.) Here is an illustration of its use within CG, or the **conjugate gradient** solver: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "it =  0  err =  0.05522230048185231\n",
      "it =  1  err =  0.09382744174195908\n",
      "it =  2  err =  0.10587866002039004\n",
      "it =  3  err =  0.09033871153608707\n",
      "it =  4  err =  0.10086392272140011\n",
      "it =  5  err =  0.08517119355086304\n",
      "it =  6  err =  0.08906575044224518\n",
      "it =  7  err =  0.07684062913929467\n",
      "it =  8  err =  0.056047175061895105\n",
      "it =  9  err =  0.04472071103380234\n",
      "it =  10  err =  0.03419741252004591\n",
      "it =  11  err =  0.02682484049645099\n",
      "it =  12  err =  0.022322983522404753\n",
      "it =  13  err =  0.020056199326844725\n",
      "it =  14  err =  0.014372559087884462\n",
      "it =  15  err =  0.010772925966275729\n",
      "it =  16  err =  0.009633149024762148\n",
      "it =  17  err =  0.0073015321836821315\n",
      "it =  18  err =  0.004811737390491996\n",
      "it =  19  err =  0.0032534005521665317\n",
      "it =  20  err =  0.001902580694234345\n",
      "it =  21  err =  0.0012857479602116576\n",
      "it =  22  err =  0.0010097069801732902\n",
      "it =  23  err =  0.0006971063184209084\n",
      "it =  24  err =  0.0004415360534270248\n",
      "it =  25  err =  0.00031313626798158353\n",
      "it =  26  err =  0.0002459345482987088\n",
      "it =  27  err =  0.00018102122685671873\n",
      "it =  28  err =  0.00012566996433136174\n",
      "it =  29  err =  9.081609378020966e-05\n",
      "it =  30  err =  5.920492010745938e-05\n",
      "it =  31  err =  4.0326988735642384e-05\n",
      "it =  32  err =  3.115865171601253e-05\n",
      "it =  33  err =  2.2328835017506844e-05\n",
      "it =  34  err =  1.63735443518241e-05\n",
      "it =  35  err =  1.1514375356105962e-05\n",
      "it =  36  err =  6.8946264410375216e-06\n",
      "it =  37  err =  4.683506893773661e-06\n",
      "it =  38  err =  3.5825089011129006e-06\n",
      "it =  39  err =  2.859816977428613e-06\n",
      "it =  40  err =  1.8589434122159337e-06\n",
      "it =  41  err =  1.2888774288482613e-06\n",
      "it =  42  err =  9.608284726114791e-07\n",
      "it =  43  err =  8.138518532407976e-07\n",
      "it =  44  err =  6.559651432751399e-07\n",
      "it =  45  err =  5.088391213652928e-07\n",
      "it =  46  err =  3.7554232878879976e-07\n",
      "it =  47  err =  2.866653641968371e-07\n",
      "it =  48  err =  2.3813379098727406e-07\n",
      "it =  49  err =  1.784776125351117e-07\n",
      "it =  50  err =  1.2713130281505307e-07\n",
      "it =  51  err =  8.879135461994969e-08\n",
      "it =  52  err =  7.189655596299818e-08\n",
      "it =  53  err =  5.3603595340687565e-08\n",
      "it =  54  err =  3.682678853570715e-08\n",
      "it =  55  err =  2.5543483270413777e-08\n",
      "it =  56  err =  1.7986625269364435e-08\n",
      "it =  57  err =  1.2658159613085641e-08\n",
      "it =  58  err =  8.74794105831229e-09\n",
      "it =  59  err =  6.004571771303868e-09\n",
      "it =  60  err =  4.240187119793316e-09\n",
      "it =  61  err =  2.9124607839179425e-09\n",
      "it =  62  err =  1.95454042937834e-09\n",
      "it =  63  err =  1.28442596904015e-09\n",
      "it =  64  err =  8.403738519702617e-10\n",
      "it =  65  err =  5.901727212281454e-10\n",
      "it =  66  err =  4.443364083343974e-10\n",
      "it =  67  err =  3.1599407725917446e-10\n",
      "it =  68  err =  2.0604768404034608e-10\n",
      "it =  69  err =  1.2843022850328627e-10\n",
      "it =  70  err =  9.063493172245921e-11\n",
      "it =  71  err =  6.766709792706894e-11\n",
      "it =  72  err =  4.746461848492383e-11\n",
      "it =  73  err =  2.7809057936360353e-11\n",
      "it =  74  err =  1.709795688455029e-11\n",
      "it =  75  err =  1.176897616349203e-11\n",
      "it =  76  err =  8.51804968461161e-12\n",
      "it =  77  err =  6.277337591794344e-12\n",
      "it =  78  err =  3.9132090932737885e-12\n",
      "it =  79  err =  2.2135436745856382e-12\n",
      "it =  80  err =  1.3679858558673982e-12\n",
      "it =  81  err =  9.755276306296867e-13\n",
      "it =  82  err =  7.137682424651503e-13\n",
      "it =  83  err =  4.4876451218901374e-13\n",
      "it =  84  err =  2.8581994514363945e-13\n",
      "it =  85  err =  1.8635801952125264e-13\n",
      "it =  86  err =  1.3655897098879978e-13\n",
      "it =  87  err =  9.36583145099722e-14\n",
      "it =  88  err =  6.003735062164544e-14\n"
     ]
    }
   ],
   "source": [
    "solvers.CG(mat=a.mat, pre=preJpoint, rhs=f.vec, sol=gfu.vec)\n",
    "Draw(gfu)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##  Gauss-Seidel smoothing\n",
    "\n",
    "The *same* point Jacobi smoother object can also used to perform **point Gauss-Seidel** smoothing. One step of the classical Gauss-Seidel iteration is realized by the method `preJpoint.Smooth()`. Its well known that this iteration converges for matrices like $A$. Below we show how to use it as a linear iterative solver. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "it# 0 , res = 0.08192679624402945\n",
      "it# 1 , res = 0.07726762555864525\n",
      "it# 2 , res = 0.07377040495808554\n",
      "it# 3 , res = 0.07062220895606912\n",
      "it# 4 , res = 0.06777582421551728\n",
      "it# 5 , res = 0.06516770753241384\n",
      "it# 6 , res = 0.06275317789329397\n",
      "it# 7 , res = 0.06050089690216452\n",
      "it# 8 , res = 0.058387526692409986\n",
      "it# 9 , res = 0.056394947347621995\n",
      "it# 10 , res = 0.0545086993995133\n",
      "it# 11 , res = 0.05271702122673418\n",
      "it# 12 , res = 0.05101020777328185\n",
      "it# 13 , res = 0.04938015983450581\n",
      "it# 14 , res = 0.04782005439140267\n",
      "it# 15 , res = 0.04632409596387132\n",
      "it# 16 , res = 0.04488732463071546\n",
      "it# 17 , res = 0.04350546518198168\n",
      "it# 18 , res = 0.042174807039847335\n",
      "it# 19 , res = 0.04089210774233109\n",
      "it# 20 , res = 0.0396545147959703\n",
      "it# 21 , res = 0.03845950204060288\n",
      "it# 22 , res = 0.03730481759400025\n",
      "it# 23 , res = 0.0361884411063836\n",
      "it# 24 , res = 0.03510854854333058\n",
      "it# 25 , res = 0.03406348308437501\n",
      "it# 26 , res = 0.033051731008142335\n",
      "it# 27 , res = 0.032071901655931344\n",
      "it# 28 , res = 0.03112271073987291\n",
      "it# 29 , res = 0.030202966400259202\n",
      "it# 30 , res = 0.029311557527407415\n",
      "it# 31 , res = 0.028447443952510115\n",
      "it# 32 , res = 0.02760964818388455\n",
      "it# 33 , res = 0.02679724842336503\n",
      "it# 34 , res = 0.026009372645006548\n",
      "it# 35 , res = 0.025245193556920734\n",
      "it# 36 , res = 0.024503924298637167\n",
      "it# 37 , res = 0.02378481475221903\n",
      "it# 38 , res = 0.023087148366535692\n",
      "it# 39 , res = 0.0224102394114777\n",
      "it# 40 , res = 0.0217534305931838\n",
      "it# 41 , res = 0.021116090973105302\n",
      "it# 42 , res = 0.020497614143416644\n",
      "it# 43 , res = 0.0198974166192704\n",
      "it# 44 , res = 0.01931493641498571\n",
      "it# 45 , res = 0.018749631776714237\n",
      "it# 46 , res = 0.018200980048633027\n",
      "it# 47 , res = 0.017668476653454283\n",
      "it# 48 , res = 0.01715163417113754\n",
      "it# 49 , res = 0.016649981502262958\n",
      "it# 50 , res = 0.016163063104665278\n",
      "it# 51 , res = 0.01569043829370675\n",
      "it# 52 , res = 0.015231680598052692\n",
      "it# 53 , res = 0.014786377164053767\n",
      "it# 54 , res = 0.014354128202872912\n",
      "it# 55 , res = 0.013934546475367345\n",
      "it# 56 , res = 0.013527256810457532\n",
      "it# 57 , res = 0.013131895653335023\n",
      "it# 58 , res = 0.012748110640368406\n",
      "it# 59 , res = 0.012375560198007516\n",
      "it# 60 , res = 0.012013913163348851\n",
      "it# 61 , res = 0.011662848424338797\n",
      "it# 62 , res = 0.011322054577854098\n",
      "it# 63 , res = 0.010991229604123473\n",
      "it# 64 , res = 0.01067008055614577\n",
      "it# 65 , res = 0.010358323262924123\n",
      "it# 66 , res = 0.010055682045473519\n",
      "it# 67 , res = 0.009761889444681583\n",
      "it# 68 , res = 0.009476685960204578\n",
      "it# 69 , res = 0.009199819799668862\n",
      "it# 70 , res = 0.008931046637527158\n",
      "it# 71 , res = 0.008670129382983446\n",
      "it# 72 , res = 0.00841683795645999\n",
      "it# 73 , res = 0.008170949074129685\n",
      "it# 74 , res = 0.007932246040080845\n",
      "it# 75 , res = 0.007700518545721331\n",
      "it# 76 , res = 0.007475562476061317\n",
      "it# 77 , res = 0.0072571797225452084\n",
      "it# 78 , res = 0.0070451780021285534\n",
      "it# 79 , res = 0.006839370682320229\n",
      "it# 80 , res = 0.006639576611929876\n",
      "it# 81 , res = 0.006445619957279532\n",
      "it# 82 , res = 0.006257330043655511\n",
      "it# 83 , res = 0.006074541201789725\n",
      "it# 84 , res = 0.005897092619174911\n",
      "it# 85 , res = 0.005724828196028781\n",
      "it# 86 , res = 0.005557596405733547\n",
      "it# 87 , res = 0.005395250159586561\n",
      "it# 88 , res = 0.005237646675707954\n",
      "it# 89 , res = 0.005084647351957455\n",
      "it# 90 , res = 0.004936117642721644\n",
      "it# 91 , res = 0.004791926939438961\n",
      "it# 92 , res = 0.004651948454736169\n",
      "it# 93 , res = 0.004516059110055665\n",
      "it# 94 , res = 0.004384139426658663\n",
      "it# 95 , res = 0.004256073419894378\n",
      "it# 96 , res = 0.004131748496628999\n",
      "it# 97 , res = 0.004011055355734409\n",
      "it# 98 , res = 0.003893887891538068\n",
      "it# 99 , res = 0.0037801431001420736\n",
      "it# 100 , res = 0.003669720988521344\n",
      "it# 101 , res = 0.003562524486313286\n",
      "it# 102 , res = 0.003458459360217854\n",
      "it# 103 , res = 0.0033574341309254496\n",
      "it# 104 , res = 0.003259359992497481\n",
      "it# 105 , res = 0.0031641507341227644\n",
      "it# 106 , res = 0.0030717226641792865\n",
      "it# 107 , res = 0.002981994536530998\n",
      "it# 108 , res = 0.0028948874789913454\n",
      "it# 109 , res = 0.002810324923890657\n",
      "it# 110 , res = 0.0027282325406816797\n",
      "it# 111 , res = 0.002648538170524317\n",
      "it# 112 , res = 0.0025711717627896324\n",
      "it# 113 , res = 0.002496065313426207\n",
      "it# 114 , res = 0.0024231528051331455\n",
      "it# 115 , res = 0.0023523701492869746\n",
      "it# 116 , res = 0.0022836551295689427\n",
      "it# 117 , res = 0.0022169473472448743\n",
      "it# 118 , res = 0.0021521881680458556\n",
      "it# 119 , res = 0.002089320670604589\n",
      "it# 120 , res = 0.002028289596400842\n",
      "it# 121 , res = 0.0019690413011714345\n",
      "it# 122 , res = 0.0019115237077420857\n",
      "it# 123 , res = 0.0018556862602383012\n",
      "it# 124 , res = 0.0018014798796365335\n",
      "it# 125 , res = 0.001748856920614219\n",
      "it# 126 , res = 0.0016977711296619743\n",
      "it# 127 , res = 0.0016481776044206479\n",
      "it# 128 , res = 0.001600032754206714\n",
      "it# 129 , res = 0.0015532942616918012\n",
      "it# 130 , res = 0.001507921045702434\n",
      "it# 131 , res = 0.0014638732251068592\n",
      "it# 132 , res = 0.0014211120837572047\n",
      "it# 133 , res = 0.0013796000364567705\n",
      "it# 134 , res = 0.001339300595921035\n",
      "it# 135 , res = 0.0013001783407047141\n",
      "it# 136 , res = 0.0012621988840657815\n",
      "it# 137 , res = 0.001225328843739243\n",
      "it# 138 , res = 0.0011895358125941567\n",
      "it# 139 , res = 0.0011547883301479946\n",
      "it# 140 , res = 0.0011210558549129947\n",
      "it# 141 , res = 0.0010883087375511042\n",
      "it# 142 , res = 0.0010565181948121436\n",
      "it# 143 , res = 0.0010256562842346165\n",
      "it# 144 , res = 0.0009956958795847124\n",
      "it# 145 , res = 0.0009666106470129384\n",
      "it# 146 , res = 0.0009383750219080032\n",
      "it# 147 , res = 0.0009109641864258245\n",
      "it# 148 , res = 0.0008843540476760172\n",
      "it# 149 , res = 0.0008585212165446786\n",
      "it# 150 , res = 0.0008334429871366016\n",
      "it# 151 , res = 0.000809097316817568\n",
      "it# 152 , res = 0.0007854628068396603\n",
      "it# 153 , res = 0.0007625186835328571\n",
      "it# 154 , res = 0.000740244780045596\n",
      "it# 155 , res = 0.0007186215186194956\n",
      "it# 156 , res = 0.0006976298933805504\n",
      "it# 157 , res = 0.0006772514536345653\n",
      "it# 158 , res = 0.0006574682876495266\n",
      "it# 159 , res = 0.0006382630069117587\n",
      "it# 160 , res = 0.0006196187308428548\n",
      "it# 161 , res = 0.0006015190719613711\n",
      "it# 162 , res = 0.0005839481214803695\n",
      "it# 163 , res = 0.0005668904353229085\n",
      "it# 164 , res = 0.0005503310205485652\n",
      "it# 165 , res = 0.0005342553221744136\n",
      "it# 166 , res = 0.0005186492103826856\n",
      "it# 167 , res = 0.0005034989681007053\n",
      "it# 168 , res = 0.000488791278944353\n",
      "it# 169 , res = 0.0004745132155138423\n",
      "it# 170 , res = 0.0004606522280309772\n",
      "it# 171 , res = 0.00044719613330827107\n",
      "it# 172 , res = 0.0004341331040406897\n",
      "it# 173 , res = 0.0004214516584100991\n",
      "it# 174 , res = 0.0004091406499930594\n",
      "it# 175 , res = 0.0003971892579639203\n",
      "it# 176 , res = 0.0003855869775833213\n",
      "it# 177 , res = 0.0003743236109654287\n",
      "it# 178 , res = 0.000363389258114413\n",
      "it# 179 , res = 0.00035277430822268883\n",
      "it# 180 , res = 0.0003424694312235002\n",
      "it# 181 , res = 0.00033246556959025305\n",
      "it# 182 , res = 0.0003227539303752261\n",
      "it# 183 , res = 0.00031332597748125803\n",
      "it# 184 , res = 0.0003041734241585562\n",
      "it# 185 , res = 0.00029528822572135686\n",
      "it# 186 , res = 0.000286662572476714\n",
      "it# 187 , res = 0.0002782888828602794\n",
      "it# 188 , res = 0.00027015979677244453\n",
      "it# 189 , res = 0.00026226816910930074\n",
      "it# 190 , res = 0.00025460706348211725\n",
      "it# 191 , res = 0.0002471697461206412\n",
      "it# 192 , res = 0.00023994967995462916\n",
      "it# 193 , res = 0.000232940518867916\n",
      "it# 194 , res = 0.00022613610212036573\n",
      "it# 195 , res = 0.00021953044893315316\n",
      "it# 196 , res = 0.00021311775323138693\n",
      "it# 197 , res = 0.0002068923785416417\n",
      "it# 198 , res = 0.00020084885303698165\n",
      "it# 199 , res = 0.0001949818647277607\n",
      "it# 200 , res = 0.00018928625679283826\n",
      "it# 201 , res = 0.00018375702304660783\n",
      "it# 202 , res = 0.00017838930353898983\n",
      "it# 203 , res = 0.0001731783802837913\n",
      "it# 204 , res = 0.0001681196731116986\n",
      "it# 205 , res = 0.00016320873564443518\n",
      "it# 206 , res = 0.00015844125138680747\n",
      "it# 207 , res = 0.00015381302993245984\n",
      "it# 208 , res = 0.00014932000328108505\n",
      "it# 209 , res = 0.00014495822226236616\n",
      "it# 210 , res = 0.00014072385306511352\n",
      "it# 211 , res = 0.0001366131738674067\n",
      "it# 212 , res = 0.00013262257156536697\n",
      "it# 213 , res = 0.00012874853859732883\n",
      "it# 214 , res = 0.0001249876698608756\n",
      "it# 215 , res = 0.00012133665971992045\n",
      "it# 216 , res = 0.0001177922990992615\n",
      "it# 217 , res = 0.00011435147266356294\n",
      "it# 218 , res = 0.00011101115607990808\n",
      "it# 219 , res = 0.00010776841335868004\n",
      "it# 220 , res = 0.00010462039427341043\n",
      "it# 221 , res = 0.00010156433185559654\n",
      "it# 222 , res = 9.859753996242925e-05\n",
      "it# 223 , res = 9.571741091609476e-05\n",
      "it# 224 , res = 9.292141321128581e-05\n",
      "it# 225 , res = 9.020708929067876e-05\n",
      "it# 226 , res = 8.757205338429242e-05\n",
      "it# 227 , res = 8.501398941303418e-05\n",
      "it# 228 , res = 8.253064895263055e-05\n",
      "it# 229 , res = 8.011984925731201e-05\n",
      "it# 230 , res = 7.777947134147551e-05\n",
      "it# 231 , res = 7.550745811731941e-05\n",
      "it# 232 , res = 7.3301812586311e-05\n",
      "it# 233 , res = 7.11605960842311e-05\n",
      "it# 234 , res = 6.908192657717135e-05\n",
      "it# 235 , res = 6.706397700723422e-05\n",
      "it# 236 , res = 6.510497368653672e-05\n",
      "it# 237 , res = 6.320319473849986e-05\n",
      "it# 238 , res = 6.135696858402421e-05\n",
      "it# 239 , res = 5.9564672472634836e-05\n",
      "it# 240 , res = 5.782473105574829e-05\n",
      "it# 241 , res = 5.613561500249084e-05\n",
      "it# 242 , res = 5.449583965495543e-05\n",
      "it# 243 , res = 5.290396372371141e-05\n",
      "it# 244 , res = 5.135858802061855e-05\n",
      "it# 245 , res = 4.9858354229300375e-05\n",
      "it# 246 , res = 4.840194371106144e-05\n",
      "it# 247 , res = 4.698807634598913e-05\n",
      "it# 248 , res = 4.561550940747264e-05\n",
      "it# 249 , res = 4.428303647042051e-05\n",
      "it# 250 , res = 4.2989486350428265e-05\n",
      "it# 251 , res = 4.173372207456912e-05\n",
      "it# 252 , res = 4.0514639881982265e-05\n",
      "it# 253 , res = 3.933116825356327e-05\n",
      "it# 254 , res = 3.818226697063176e-05\n",
      "it# 255 , res = 3.706692620003855e-05\n",
      "it# 256 , res = 3.5984165606928974e-05\n",
      "it# 257 , res = 3.493303349294157e-05\n",
      "it# 258 , res = 3.391260595975471e-05\n",
      "it# 259 , res = 3.2921986097031015e-05\n",
      "it# 260 , res = 3.196030319400556e-05\n",
      "it# 261 , res = 3.102671197426804e-05\n",
      "it# 262 , res = 3.012039185265958e-05\n",
      "it# 263 , res = 2.924054621418891e-05\n",
      "it# 264 , res = 2.83864017136738e-05\n",
      "it# 265 , res = 2.7557207596285808e-05\n",
      "it# 266 , res = 2.6752235037222646e-05\n",
      "it# 267 , res = 2.597077650136112e-05\n",
      "it# 268 , res = 2.5212145121445153e-05\n",
      "it# 269 , res = 2.4475674094500226e-05\n",
      "it# 270 , res = 2.3760716095100676e-05\n",
      "it# 271 , res = 2.3066642707015354e-05\n",
      "it# 272 , res = 2.2392843870726235e-05\n",
      "it# 273 , res = 2.173872734695419e-05\n",
      "it# 274 , res = 2.1103718196545164e-05\n",
      "it# 275 , res = 2.048725827450787e-05\n",
      "it# 276 , res = 1.9888805740331945e-05\n",
      "it# 277 , res = 1.9307834580586206e-05\n",
      "it# 278 , res = 1.8743834147660917e-05\n",
      "it# 279 , res = 1.8196308710329865e-05\n",
      "it# 280 , res = 1.7664777018027522e-05\n",
      "it# 281 , res = 1.7148771878047257e-05\n",
      "it# 282 , res = 1.6647839744876715e-05\n",
      "it# 283 , res = 1.6161540321541195e-05\n",
      "it# 284 , res = 1.568944617237857e-05\n",
      "it# 285 , res = 1.5231142347775923e-05\n",
      "it# 286 , res = 1.4786226019085057e-05\n",
      "it# 287 , res = 1.4354306124536763e-05\n",
      "it# 288 , res = 1.3935003025872397e-05\n",
      "it# 289 , res = 1.3527948174279261e-05\n",
      "it# 290 , res = 1.3132783786710763e-05\n",
      "it# 291 , res = 1.2749162531343342e-05\n",
      "it# 292 , res = 1.2376747222128253e-05\n",
      "it# 293 , res = 1.2015210522471506e-05\n",
      "it# 294 , res = 1.1664234657852027e-05\n",
      "it# 295 , res = 1.1323511136057308e-05\n",
      "it# 296 , res = 1.0992740476262984e-05\n",
      "it# 297 , res = 1.0671631945903052e-05\n",
      "it# 298 , res = 1.035990330474014e-05\n",
      "it# 299 , res = 1.005728055716698e-05\n",
      "it# 300 , res = 9.763497711510873e-06\n",
      "it# 301 , res = 9.478296545492734e-06\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "it# 302 , res = 9.20142638007839e-06\n",
      "it# 303 , res = 8.93264385864601e-06\n",
      "it# 304 , res = 8.671712733319411e-06\n",
      "it# 305 , res = 8.418403657416411e-06\n",
      "it# 306 , res = 8.172493983483656e-06\n",
      "it# 307 , res = 7.933767567766574e-06\n",
      "it# 308 , res = 7.702014580399288e-06\n",
      "it# 309 , res = 7.477031320963506e-06\n",
      "it# 310 , res = 7.25862003899917e-06\n",
      "it# 311 , res = 7.046588760873389e-06\n",
      "it# 312 , res = 6.84075112049103e-06\n",
      "it# 313 , res = 6.640926195685272e-06\n",
      "it# 314 , res = 6.446938349203836e-06\n",
      "it# 315 , res = 6.25861707457476e-06\n",
      "it# 316 , res = 6.075796845556455e-06\n",
      "it# 317 , res = 5.898316971400937e-06\n",
      "it# 318 , res = 5.72602145524725e-06\n",
      "it# 319 , res = 5.5587588569724775e-06\n",
      "it# 320 , res = 5.396382160124041e-06\n",
      "it# 321 , res = 5.2387486429536435e-06\n",
      "it# 322 , res = 5.085719752475799e-06\n",
      "it# 323 , res = 4.9371609832884464e-06\n",
      "it# 324 , res = 4.792941758633401e-06\n",
      "it# 325 , res = 4.652935316508108e-06\n",
      "it# 326 , res = 4.5170185971830715e-06\n",
      "it# 327 , res = 4.385072136057004e-06\n",
      "it# 328 , res = 4.256979958117818e-06\n",
      "it# 329 , res = 4.132629475959835e-06\n",
      "it# 330 , res = 4.011911391118441e-06\n",
      "it# 331 , res = 3.894719597675873e-06\n",
      "it# 332 , res = 3.7809510893390857e-06\n",
      "it# 333 , res = 3.670505868684667e-06\n",
      "it# 334 , res = 3.56328685916848e-06\n",
      "it# 335 , res = 3.459199820193475e-06\n",
      "it# 336 , res = 3.3581532637970477e-06\n",
      "it# 337 , res = 3.2600583745142916e-06\n",
      "it# 338 , res = 3.1648289314628957e-06\n",
      "it# 339 , res = 3.072381232041618e-06\n",
      "it# 340 , res = 2.9826340189667907e-06\n",
      "it# 341 , res = 2.8955084082790716e-06\n",
      "it# 342 , res = 2.8109278206005745e-06\n",
      "it# 343 , res = 2.7288179132925763e-06\n",
      "it# 344 , res = 2.6491065154456158e-06\n",
      "it# 345 , res = 2.571723564207027e-06\n",
      "it# 346 , res = 2.4966010434604014e-06\n",
      "it# 347 , res = 2.423672923871059e-06\n",
      "it# 348 , res = 2.3528751048535465e-06\n",
      "it# 349 , res = 2.2841453582609813e-06\n",
      "it# 350 , res = 2.2174232737614964e-06\n",
      "it# 351 , res = 2.1526502055997703e-06\n",
      "it# 352 , res = 2.0897692210266675e-06\n",
      "it# 353 , res = 2.0287250505702893e-06\n",
      "it# 354 , res = 1.969464039117375e-06\n",
      "it# 355 , res = 1.911934098781118e-06\n",
      "it# 356 , res = 1.8560846634527472e-06\n",
      "it# 357 , res = 1.8018666438437067e-06\n",
      "it# 358 , res = 1.7492323847821698e-06\n",
      "it# 359 , res = 1.6981356230441588e-06\n",
      "it# 360 , res = 1.6485314469923352e-06\n",
      "it# 361 , res = 1.600376256649062e-06\n",
      "it# 362 , res = 1.5536277256803486e-06\n",
      "it# 363 , res = 1.508244764407703e-06\n",
      "it# 364 , res = 1.4641874830384717e-06\n",
      "it# 365 , res = 1.421417157300543e-06\n",
      "it# 366 , res = 1.379896194015483e-06\n",
      "it# 367 , res = 1.3395880979599996e-06\n",
      "it# 368 , res = 1.3004574401699742e-06\n",
      "it# 369 , res = 1.2624698266073e-06\n",
      "it# 370 , res = 1.2255918677765647e-06\n",
      "it# 371 , res = 1.1897911495989067e-06\n",
      "it# 372 , res = 1.1550362048340911e-06\n",
      "it# 373 , res = 1.1212964855609552e-06\n",
      "it# 374 , res = 1.0885423358983803e-06\n",
      "it# 375 , res = 1.0567449665072173e-06\n",
      "it# 376 , res = 1.025876428871224e-06\n",
      "it# 377 , res = 9.959095908646414e-07\n",
      "it# 378 , res = 9.668181131464157e-07\n",
      "it# 379 , res = 9.385764254222216e-07\n",
      "it# 380 , res = 9.111597045982461e-07\n",
      "it# 381 , res = 8.845438526144883e-07\n",
      "it# 382 , res = 8.587054752484367e-07\n",
      "it# 383 , res = 8.336218617996018e-07\n",
      "it# 384 , res = 8.092709647933546e-07\n",
      "it# 385 , res = 7.856313809619995e-07\n",
      "it# 386 , res = 7.62682332102173e-07\n",
      "it# 387 , res = 7.404036470599436e-07\n",
      "it# 388 , res = 7.187757438579928e-07\n",
      "it# 389 , res = 6.977796125536263e-07\n",
      "it# 390 , res = 6.773967983881534e-07\n",
      "it# 391 , res = 6.576093858627596e-07\n",
      "it# 392 , res = 6.38399982732111e-07\n",
      "it# 393 , res = 6.197517047638082e-07\n",
      "it# 394 , res = 6.01648160936633e-07\n",
      "it# 395 , res = 5.840734389872142e-07\n",
      "it# 396 , res = 5.670120916469611e-07\n",
      "it# 397 , res = 5.50449122686387e-07\n",
      "it# 398 , res = 5.343699739852417e-07\n",
      "it# 399 , res = 5.187605126956174e-07\n",
      "it# 400 , res = 5.036070187690837e-07\n",
      "it# 401 , res = 4.888961729790659e-07\n",
      "it# 402 , res = 4.7461504513396437e-07\n",
      "it# 403 , res = 4.6075108286727444e-07\n",
      "it# 404 , res = 4.472921002565832e-07\n",
      "it# 405 , res = 4.3422626753263947e-07\n",
      "it# 406 , res = 4.21542100204526e-07\n",
      "it# 407 , res = 4.092284497445619e-07\n",
      "it# 408 , res = 3.9727449282075296e-07\n",
      "it# 409 , res = 3.8566972244087746e-07\n",
      "it# 410 , res = 3.7440393860630063e-07\n",
      "it# 411 , res = 3.6346723919561035e-07\n",
      "it# 412 , res = 3.5285001131485737e-07\n",
      "it# 413 , res = 3.4254292287366425e-07\n",
      "it# 414 , res = 3.325369143416729e-07\n",
      "it# 415 , res = 3.2282319100270845e-07\n",
      "it# 416 , res = 3.13393214867448e-07\n",
      "it# 417 , res = 3.0423869744928535e-07\n",
      "it# 418 , res = 2.953515921734483e-07\n",
      "it# 419 , res = 2.8672408788377657e-07\n",
      "it# 420 , res = 2.7834860134073113e-07\n",
      "it# 421 , res = 2.702177706588409e-07\n",
      "it# 422 , res = 2.6232444952714855e-07\n",
      "it# 423 , res = 2.5466169985204204e-07\n",
      "it# 424 , res = 2.4722278645321504e-07\n",
      "it# 425 , res = 2.400011709010443e-07\n",
      "it# 426 , res = 2.329905057260869e-07\n",
      "it# 427 , res = 2.2618462877411483e-07\n",
      "it# 428 , res = 2.1957755802810812e-07\n",
      "it# 429 , res = 2.131634862099809e-07\n",
      "it# 430 , res = 2.069367755911583e-07\n",
      "it# 431 , res = 2.0089195322001486e-07\n",
      "it# 432 , res = 1.9502370583299227e-07\n",
      "it# 433 , res = 1.8932687573653563e-07\n",
      "it# 434 , res = 1.8379645543772082e-07\n",
      "it# 435 , res = 1.7842758409160176e-07\n",
      "it# 436 , res = 1.7321554260859983e-07\n",
      "it# 437 , res = 1.6815574985223712e-07\n",
      "it# 438 , res = 1.6324375852972194e-07\n",
      "it# 439 , res = 1.5847525120331187e-07\n",
      "it# 440 , res = 1.538460365826503e-07\n",
      "it# 441 , res = 1.4935204577373734e-07\n",
      "it# 442 , res = 1.449893287671978e-07\n",
      "it# 443 , res = 1.4075405091746439e-07\n",
      "it# 444 , res = 1.366424896713164e-07\n",
      "it# 445 , res = 1.3265103096566273e-07\n",
      "it# 446 , res = 1.2877616666977125e-07\n",
      "it# 447 , res = 1.250144908844406e-07\n",
      "it# 448 , res = 1.2136269731172844e-07\n",
      "it# 449 , res = 1.1781757605355208e-07\n",
      "it# 450 , res = 1.1437601128620873e-07\n",
      "it# 451 , res = 1.1103497789963403e-07\n",
      "it# 452 , res = 1.0779153926650003e-07\n",
      "it# 453 , res = 1.0464284469397545e-07\n",
      "it# 454 , res = 1.0158612644480159e-07\n",
      "it# 455 , res = 9.861869787720146e-08\n",
      "it# 456 , res = 9.573795077326162e-08\n",
      "it# 457 , res = 9.294135313771748e-08\n",
      "it# 458 , res = 9.022644680141832e-08\n",
      "it# 459 , res = 8.759084540422157e-08\n",
      "it# 460 , res = 8.503223247641311e-08\n",
      "it# 461 , res = 8.254835915869222e-08\n",
      "it# 462 , res = 8.01370421416458e-08\n",
      "it# 463 , res = 7.77961620073764e-08\n",
      "it# 464 , res = 7.552366120809265e-08\n",
      "it# 465 , res = 7.331754236129858e-08\n",
      "it# 466 , res = 7.117586639524648e-08\n",
      "it# 467 , res = 6.909675081918946e-08\n",
      "it# 468 , res = 6.707836818986148e-08\n",
      "it# 469 , res = 6.511894450425933e-08\n",
      "it# 470 , res = 6.321675746324933e-08\n",
      "it# 471 , res = 6.137013513206768e-08\n",
      "it# 472 , res = 5.957745439495801e-08\n",
      "it# 473 , res = 5.7837139561921364e-08\n",
      "it# 474 , res = 5.614766108888818e-08\n",
      "it# 475 , res = 5.4507533852471614e-08\n",
      "it# 476 , res = 5.291531629186205e-08\n",
      "it# 477 , res = 5.1369608995322195e-08\n",
      "it# 478 , res = 4.986905328823826e-08\n",
      "it# 479 , res = 4.841233020610855e-08\n",
      "it# 480 , res = 4.699815949500824e-08\n",
      "it# 481 , res = 4.56252979541951e-08\n",
      "it# 482 , res = 4.4292539149818404e-08\n",
      "it# 483 , res = 4.299871142157872e-08\n",
      "it# 484 , res = 4.174267767994465e-08\n",
      "it# 485 , res = 4.052333390190662e-08\n",
      "it# 486 , res = 3.9339608322299644e-08\n",
      "it# 487 , res = 3.819046044416984e-08\n",
      "it# 488 , res = 3.7074880325813256e-08\n",
      "it# 489 , res = 3.5991887391963725e-08\n",
      "it# 490 , res = 3.4940529730331197e-08\n",
      "it# 491 , res = 3.391988323546586e-08\n",
      "it# 492 , res = 3.292905081126291e-08\n",
      "it# 493 , res = 3.196716148699039e-08\n",
      "it# 494 , res = 3.1033369958910765e-08\n",
      "it# 495 , res = 3.012685531488535e-08\n",
      "it# 496 , res = 2.9246820885573157e-08\n",
      "it# 497 , res = 2.8392493149548758e-08\n",
      "it# 498 , res = 2.7563121083519314e-08\n",
      "it# 499 , res = 2.6757975813604662e-08\n"
     ]
    }
   ],
   "source": [
    "gfu.vec[:] = 0\n",
    "res = f.vec.CreateVector()              # residual \n",
    "projres = f.vec.CreateVector()          # residual projected to freedofs\n",
    "proj = Projector(fes.FreeDofs(), True)\n",
    "\n",
    "for i in range(500):\n",
    "    preJpoint.Smooth(gfu.vec, f.vec)    # one step of point Gauss-Seidel\n",
    "    res.data = f.vec - a.mat*gfu.vec      \n",
    "    projres.data = proj * res\n",
    "    print (\"it#\", i, \", res =\", Norm(projres))\n",
    "Draw (gfu)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Implement a forward-backward GS preconditioner"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The *same* point Jacobi smoother object is also able to perform a Gauss-Seidel iteration after reversing the ordering of the points, i.e., a **backward** Gauss-Seidel sweep. One can combine the forward and backward sweeps to construct a symmetric preconditioner, often called the **symmetrized Gauss-Seidel preconditioner**. This offers a good illustration of how to construct NGSolve preconditioners from within python. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "class SymmetricGS(BaseMatrix):\n",
    "    def __init__ (self, smoother):\n",
    "        super(SymmetricGS, self).__init__()\n",
    "        self.smoother = smoother\n",
    "    def Mult (self, x, y):\n",
    "        y[:] = 0.0\n",
    "        self.smoother.Smooth(y, x)\n",
    "        self.smoother.SmoothBack(y,x)\n",
    "    def Height (self):\n",
    "        return self.smoother.height\n",
    "    def Width (self):\n",
    "        return self.smoother.height"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "it =  0  err =  0.09421429683611413\n",
      "it =  1  err =  0.12242092674801901\n",
      "it =  2  err =  0.08627911974342342\n",
      "it =  3  err =  0.054294636877542546\n",
      "it =  4  err =  0.026834953641369612\n",
      "it =  5  err =  0.014442946534311274\n",
      "it =  6  err =  0.006828162088958496\n",
      "it =  7  err =  0.002947538932807923\n",
      "it =  8  err =  0.0009884874866274685\n",
      "it =  9  err =  0.00037122112200935675\n",
      "it =  10  err =  0.00011770383684169319\n",
      "it =  11  err =  4.7294009691925314e-05\n",
      "it =  12  err =  1.795956404345207e-05\n",
      "it =  13  err =  7.73981616121209e-06\n",
      "it =  14  err =  3.688537792536532e-06\n",
      "it =  15  err =  1.7853513785615997e-06\n",
      "it =  16  err =  6.793337817618327e-07\n",
      "it =  17  err =  3.267271476258251e-07\n",
      "it =  18  err =  1.0509766947753396e-07\n",
      "it =  19  err =  5.2632976611271436e-08\n",
      "it =  20  err =  1.7518834053967413e-08\n",
      "it =  21  err =  7.1749421044469985e-09\n",
      "it =  22  err =  2.5718606236788454e-09\n",
      "it =  23  err =  7.674623733376304e-10\n",
      "it =  24  err =  3.2317029440816587e-10\n",
      "it =  25  err =  7.901568725512893e-11\n",
      "it =  26  err =  2.830256119137346e-11\n",
      "it =  27  err =  1.0594770382019746e-11\n",
      "it =  28  err =  2.983199288638775e-12\n",
      "it =  29  err =  1.2174564854809225e-12\n",
      "it =  30  err =  4.3306940793653186e-13\n",
      "it =  31  err =  1.6903295709170732e-13\n"
     ]
    }
   ],
   "source": [
    "preGS = SymmetricGS(preJpoint)\n",
    "solvers.CG(mat=a.mat, pre=preGS, rhs=f.vec, sol=gfu.vec)\n",
    "Draw (gfu)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "20.165532778155335"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lams = EigenValues_Preconditioner(mat=a.mat, pre=preGS)\n",
    "max(lams)/min(lams)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that the condition number now is better than that of the system preconditioned by point Jacobi."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## A Block Jacobi preconditioner\n",
    "\n",
    "The point Jacobi preconditioner is based on inverses of 1 x 1 diagonal blocks.  Condition numbers can be improved by using larger blocks. It is possible to group dofs into blocks within python and construct an NGSolve preconditioner based on the blocks.\n",
    "\n",
    "Here is an example that constructs vertex-based blocks."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[{866, 158, 159}, {13, 206, 207, 48, 145, 144, 882, 142, 143, 210, 211, 884}, {256, 257, 2, 901, 146, 147, 148, 21, 22, 149}, {65, 314, 919, 315, 917, 310, 150, 151, 154, 155, 311, 30}, {160, 161, 866, 867, 164, 165, 935, 40, 361, 360, 158, 159}, {869, 160, 161, 867, 164, 165, 868, 166, 40, 167, 41, 170, 171, 362, 363}, {868, 166, 167, 870, 41, 170, 171, 42, 172, 173, 871, 176, 177, 368, 369}, {375, 870, 872, 873, 42, 43, 172, 173, 176, 177, 178, 179, 182, 183, 374}, {380, 872, 874, 43, 44, 875, 178, 179, 381, 182, 183, 184, 185, 188, 189}, {194, 195, 386, 387, 874, 44, 876, 45, 877, 184, 185, 188, 189, 190, 191}, {194, 195, 196, 197, 200, 201, 392, 393, 876, 45, 46, 878, 879, 190, 191}, {196, 197, 200, 201, 202, 203, 204, 205, 398, 399, 878, 46, 47, 880, 881}, {202, 203, 204, 205, 206, 207, 144, 145, 404, 405, 47, 880, 48, 882, 883}, {13, 142, 143, 144, 145, 210, 211, 14, 208, 209, 212, 213, 217, 216, 410, 411, 48, 49, 884, 885, 886}, {13, 14, 15, 208, 209, 212, 213, 214, 215, 216, 217, 218, 219, 222, 223, 414, 415, 49, 50, 885, 887, 888}, {14, 15, 16, 214, 215, 218, 219, 220, 221, 222, 223, 224, 225, 228, 229, 420, 421, 50, 51, 887, 889, 890}, {15, 16, 17, 220, 221, 224, 225, 226, 227, 228, 229, 230, 231, 234, 235, 426, 427, 51, 52, 889, 891, 892}, {16, 17, 18, 226, 227, 230, 231, 232, 233, 234, 235, 236, 237, 240, 241, 432, 433, 52, 53, 891, 893, 894}, {896, 17, 18, 19, 247, 439, 232, 233, 236, 237, 238, 239, 240, 241, 242, 243, 53, 54, 246, 438, 893, 895}, {897, 898, 18, 19, 20, 55, 247, 445, 238, 239, 242, 243, 244, 245, 54, 246, 248, 249, 444, 252, 253, 895}, {897, 258, 259, 899, 450, 451, 900, 19, 20, 21, 248, 244, 245, 55, 56, 249, 250, 251, 252, 253, 254, 255}, {256, 257, 258, 259, 899, 2, 901, 262, 263, 146, 147, 20, 21, 148, 22, 149, 936, 56, 250, 251, 254, 255}, {256, 257, 2, 258, 260, 901, 261, 902, 264, 265, 259, 262, 268, 269, 263, 456, 457, 146, 147, 148, 21, 22, 149, 23, 936, 937, 56, 57}, {260, 261, 902, 903, 264, 265, 266, 267, 268, 269, 270, 271, 904, 460, 274, 275, 461, 22, 23, 24, 57, 58}, {903, 905, 266, 267, 906, 270, 271, 272, 273, 274, 275, 276, 277, 464, 23, 280, 24, 281, 25, 465, 58, 59}, {471, 905, 907, 908, 272, 273, 276, 277, 278, 279, 280, 281, 26, 283, 24, 25, 282, 286, 287, 59, 60, 470}, {476, 907, 909, 910, 278, 279, 25, 282, 27, 284, 26, 283, 285, 288, 289, 286, 287, 292, 293, 477, 60, 61}, {909, 911, 912, 26, 27, 28, 285, 284, 288, 289, 290, 291, 292, 293, 294, 295, 482, 483, 298, 299, 61, 62}, {911, 913, 914, 27, 28, 29, 290, 291, 294, 295, 296, 297, 298, 299, 300, 301, 488, 489, 304, 305, 62, 63}, {64, 913, 915, 916, 28, 29, 30, 296, 297, 300, 301, 302, 303, 304, 305, 306, 307, 308, 309, 494, 495, 63}, {64, 65, 915, 917, 150, 151, 918, 154, 155, 29, 30, 302, 303, 306, 307, 308, 309, 310, 311, 500, 501}, {320, 65, 321, 66, 919, 920, 921, 154, 155, 506, 507, 314, 315, 316, 317}, {320, 321, 66, 322, 323, 67, 326, 327, 920, 922, 923, 316, 317, 510, 511}, {322, 323, 67, 68, 326, 327, 328, 329, 516, 517, 332, 333, 922, 924, 925}, {68, 69, 328, 329, 522, 523, 332, 333, 334, 335, 338, 339, 924, 926, 927}, {69, 70, 334, 335, 528, 529, 338, 339, 340, 341, 344, 345, 926, 928, 929}, {70, 71, 340, 341, 534, 535, 344, 345, 346, 347, 350, 351, 928, 930, 931}, {71, 72, 346, 347, 540, 541, 350, 351, 352, 353, 930, 932, 933, 358, 359}, {72, 352, 353, 932, 356, 358, 359, 357, 934, 40, 361, 360, 938, 364, 365}, {160, 161, 866, 356, 357, 934, 935, 40, 361, 360, 158, 159}, {158, 159, 160, 161, 546, 547, 164, 165, 166, 935, 167, 40, 41, 934, 938, 940, 999, 72, 74, 867, 356, 869, 358, 359, 357, 361, 362, 363, 360, 364, 365, 366, 367, 372, 373}, {73, 74, 550, 367, 869, 551, 868, 164, 166, 167, 165, 41, 170, 171, 40, 362, 363, 871, 42, 172, 173, 368, 369, 939, 370, 371, 376, 377, 940, 366, 372, 373, 941}, {384, 385, 73, 377, 997, 552, 375, 553, 994, 376, 100, 870, 871, 378, 41, 42, 170, 172, 173, 171, 873, 176, 177, 368, 369, 43, 178, 179, 374, 939, 370, 371, 379}, {384, 385, 388, 389, 75, 184, 375, 994, 995, 100, 872, 873, 42, 43, 44, 875, 942, 176, 177, 178, 179, 564, 565, 182, 183, 374, 185, 378, 379, 380, 381, 382, 383}, {386, 387, 388, 389, 382, 383, 390, 391, 394, 75, 395, 76, 189, 380, 874, 43, 44, 875, 45, 877, 942, 943, 944, 562, 563, 182, 183, 184, 185, 188, 381, 190, 191}, {194, 195, 386, 387, 196, 197, 392, 393, 390, 391, 394, 395, 76, 396, 77, 397, 400, 401, 876, 45, 44, 877, 46, 879, 943, 945, 946, 570, 571, 188, 189, 190, 191}, {576, 577, 194, 195, 196, 197, 200, 201, 392, 393, 202, 203, 398, 399, 396, 77, 397, 400, 401, 78, 402, 403, 406, 407, 45, 878, 46, 879, 47, 881, 945, 947, 948}, {582, 79, 200, 201, 202, 203, 204, 205, 398, 399, 206, 207, 78, 402, 404, 405, 403, 406, 407, 408, 409, 412, 413, 46, 47, 880, 881, 48, 883, 947, 949, 583, 950}, {79, 204, 205, 206, 207, 144, 145, 13, 142, 404, 405, 143, 210, 211, 212, 213, 410, 411, 408, 409, 412, 413, 416, 417, 47, 48, 49, 882, 883, 884, 949, 886, 951}, {13, 14, 410, 411, 412, 413, 414, 415, 416, 417, 418, 419, 424, 425, 48, 49, 50, 951, 953, 954, 588, 589, 79, 208, 209, 210, 211, 212, 213, 81, 216, 217, 218, 219, 885, 886, 888}, {14, 15, 414, 415, 418, 419, 420, 421, 422, 423, 424, 425, 428, 429, 49, 50, 51, 952, 953, 955, 590, 591, 80, 81, 214, 215, 216, 217, 218, 219, 222, 223, 224, 225, 887, 888, 890}, {15, 16, 420, 421, 422, 423, 426, 427, 428, 429, 430, 431, 50, 51, 52, 434, 435, 952, 956, 957, 80, 592, 82, 593, 220, 221, 222, 223, 224, 225, 228, 229, 230, 231, 889, 890, 892}, {16, 17, 426, 427, 430, 431, 432, 433, 434, 51, 52, 53, 435, 436, 440, 441, 437, 956, 958, 959, 82, 83, 600, 601, 226, 227, 228, 229, 230, 231, 234, 235, 236, 237, 891, 892, 894}, {896, 17, 18, 432, 433, 52, 53, 54, 439, 440, 441, 438, 443, 436, 437, 446, 447, 960, 958, 442, 961, 83, 84, 606, 607, 232, 233, 234, 235, 236, 237, 240, 241, 242, 243, 893, 894}, {896, 898, 18, 19, 53, 54, 55, 439, 438, 442, 443, 444, 445, 446, 447, 960, 448, 449, 962, 452, 453, 963, 84, 85, 612, 613, 238, 239, 240, 241, 242, 243, 246, 247, 248, 249, 895}, {897, 898, 900, 19, 20, 54, 55, 56, 444, 445, 448, 449, 450, 451, 962, 452, 453, 964, 454, 455, 458, 459, 965, 85, 86, 618, 619, 244, 245, 246, 247, 248, 249, 252, 253, 254, 255}, {256, 257, 258, 259, 899, 900, 262, 263, 264, 265, 20, 21, 22, 936, 937, 55, 56, 57, 450, 451, 964, 454, 455, 456, 457, 458, 459, 966, 462, 463, 86, 250, 251, 252, 253, 254, 255}, {260, 261, 902, 262, 264, 265, 904, 263, 268, 269, 460, 461, 270, 271, 457, 458, 459, 462, 22, 23, 463, 86, 970, 466, 467, 937, 456, 56, 57, 58, 966}, {903, 904, 266, 267, 268, 269, 270, 271, 906, 274, 275, 276, 277, 23, 24, 57, 58, 59, 968, 970, 460, 461, 462, 463, 464, 465, 466, 467, 468, 469, 86, 88, 474, 475, 996, 624, 625}, {905, 906, 908, 272, 273, 274, 275, 276, 277, 280, 281, 24, 25, 283, 282, 58, 59, 60, 967, 968, 969, 464, 465, 468, 469, 470, 471, 87, 472, 473, 88, 474, 475, 478, 479, 626, 627}, {907, 908, 910, 278, 279, 280, 281, 282, 25, 26, 283, 286, 287, 288, 289, 59, 60, 61, 967, 971, 972, 470, 471, 472, 89, 473, 87, 476, 477, 478, 479, 480, 481, 484, 485, 628, 629}, {909, 910, 912, 26, 27, 284, 285, 286, 287, 288, 289, 292, 293, 294, 295, 60, 61, 62, 971, 973, 974, 89, 90, 476, 477, 480, 481, 482, 483, 484, 485, 486, 487, 490, 491, 638, 639}, {644, 645, 911, 912, 914, 27, 28, 290, 291, 292, 293, 294, 295, 298, 299, 300, 301, 61, 62, 63, 973, 975, 976, 90, 91, 482, 483, 486, 487, 488, 489, 490, 491, 492, 493, 496, 497}, {650, 651, 913, 914, 916, 28, 29, 296, 297, 298, 299, 300, 301, 304, 305, 306, 307, 62, 63, 64, 975, 977, 978, 91, 92, 488, 489, 492, 493, 494, 495, 496, 497, 498, 499, 502, 503}, {656, 657, 915, 916, 918, 29, 30, 302, 303, 304, 305, 306, 307, 308, 309, 310, 311, 63, 64, 65, 977, 979, 980, 92, 93, 494, 495, 498, 499, 500, 501, 502, 503, 504, 505, 508, 509}, {64, 65, 66, 509, 512, 513, 981, 979, 507, 917, 150, 151, 918, 919, 154, 155, 921, 311, 30, 93, 504, 508, 506, 308, 501, 310, 309, 500, 505, 314, 315, 316, 317}, {320, 321, 66, 65, 67, 322, 323, 512, 509, 513, 514, 515, 520, 521, 981, 662, 983, 920, 921, 663, 923, 984, 93, 95, 314, 508, 507, 506, 315, 316, 317, 510, 511}, {320, 321, 322, 323, 67, 66, 326, 327, 68, 516, 517, 328, 329, 518, 519, 524, 525, 520, 521, 982, 983, 664, 665, 922, 923, 514, 925, 94, 95, 515, 985, 510, 511}, {67, 68, 516, 326, 327, 328, 329, 517, 69, 332, 333, 522, 523, 334, 335, 524, 525, 526, 527, 982, 530, 531, 986, 987, 924, 925, 94, 927, 96, 666, 667, 518, 519}, {68, 69, 70, 537, 522, 523, 332, 333, 334, 335, 528, 529, 338, 339, 340, 341, 526, 527, 530, 531, 986, 532, 533, 536, 926, 927, 96, 929, 97, 988, 674, 675, 989}, {69, 70, 71, 528, 529, 338, 339, 340, 341, 534, 535, 344, 345, 346, 347, 532, 533, 536, 537, 928, 929, 97, 931, 988, 98, 990, 538, 542, 680, 681, 543, 539, 991}, {992, 993, 70, 71, 72, 538, 534, 535, 344, 345, 346, 347, 540, 541, 350, 351, 352, 353, 930, 931, 98, 933, 990, 542, 543, 544, 545, 99, 548, 549, 686, 687, 539}, {992, 547, 548, 71, 72, 74, 540, 541, 350, 351, 352, 353, 544, 545, 932, 933, 358, 359, 356, 357, 40, 938, 364, 365, 549, 998, 558, 559, 999, 99, 366, 367, 546}, {550, 551, 552, 41, 42, 939, 556, 553, 941, 554, 560, 561, 555, 557, 694, 695, 73, 74, 754, 755, 100, 997, 111, 368, 369, 370, 371, 372, 373, 379, 112, 376, 377, 378, 1019, 1020, 1022}, {546, 547, 548, 549, 550, 551, 40, 41, 940, 941, 558, 559, 560, 556, 557, 561, 692, 693, 72, 73, 74, 99, 998, 999, 362, 363, 364, 365, 366, 367, 112, 370, 371, 372, 373, 1020, 1021}, {384, 385, 388, 389, 390, 391, 1035, 1036, 43, 44, 942, 944, 562, 563, 564, 565, 566, 567, 696, 568, 569, 697, 574, 575, 708, 709, 75, 76, 995, 100, 102, 1001, 118, 380, 381, 382, 383}, {386, 387, 388, 389, 390, 391, 394, 395, 396, 397, 44, 45, 943, 944, 562, 563, 946, 566, 567, 570, 571, 572, 573, 574, 575, 698, 699, 578, 579, 75, 76, 77, 101, 102, 1000, 1001, 1002}, {392, 393, 394, 395, 396, 397, 1037, 400, 401, 402, 403, 45, 46, 945, 946, 948, 570, 571, 572, 573, 700, 701, 576, 577, 578, 579, 580, 581, 586, 587, 76, 77, 78, 101, 103, 1000, 1018}, {398, 399, 400, 401, 402, 403, 406, 407, 408, 409, 46, 47, 947, 948, 950, 576, 577, 580, 581, 582, 583, 584, 585, 586, 587, 588, 77, 78, 79, 589, 81, 598, 599, 103, 1003, 1004, 1018}, {582, 583, 584, 585, 588, 589, 78, 79, 81, 404, 405, 406, 407, 408, 409, 410, 411, 412, 413, 416, 417, 418, 419, 1003, 47, 48, 49, 949, 950, 951, 954}, {1038, 1040, 420, 421, 422, 423, 424, 425, 428, 429, 430, 431, 50, 51, 952, 955, 957, 714, 715, 590, 591, 80, 81, 592, 82, 593, 596, 597, 594, 595, 598, 599, 604, 605, 103, 106, 1008}, {1038, 414, 415, 416, 417, 418, 419, 422, 423, 424, 425, 49, 50, 953, 954, 955, 582, 583, 584, 585, 586, 587, 588, 589, 590, 79, 591, 81, 80, 78, 594, 595, 598, 599, 103, 1003, 1004}, {426, 427, 428, 429, 430, 431, 434, 51, 52, 435, 436, 437, 956, 957, 959, 718, 719, 80, 592, 82, 593, 83, 596, 597, 600, 601, 602, 603, 604, 605, 608, 609, 104, 106, 1005, 1008, 1009}, {432, 433, 434, 435, 436, 52, 53, 437, 440, 441, 442, 443, 958, 959, 961, 716, 717, 82, 83, 84, 600, 601, 602, 603, 606, 607, 608, 609, 610, 611, 614, 615, 104, 105, 1005, 1006, 1007}, {53, 54, 439, 440, 438, 441, 442, 443, 446, 447, 960, 961, 448, 449, 963, 83, 84, 85, 724, 725, 606, 607, 610, 611, 612, 613, 614, 615, 616, 105, 617, 107, 1006, 622, 623, 1010, 1011}, {1034, 54, 55, 444, 445, 446, 447, 448, 449, 962, 963, 452, 453, 454, 455, 965, 84, 85, 86, 88, 612, 613, 616, 617, 618, 619, 107, 620, 622, 623, 621, 624, 1010, 625, 1012, 634, 635}, {55, 56, 57, 58, 450, 451, 964, 452, 454, 455, 453, 965, 458, 459, 966, 456, 457, 462, 463, 970, 460, 461, 466, 85, 86, 467, 468, 469, 88, 996, 618, 619, 620, 621, 624, 625, 1012}, {642, 643, 59, 60, 1016, 967, 969, 972, 470, 471, 88, 87, 474, 472, 473, 475, 478, 479, 480, 481, 89, 744, 745, 109, 110, 633, 626, 627, 628, 629, 630, 631, 632, 1014, 1017, 636, 637}, {1034, 1041, 107, 58, 59, 968, 969, 464, 465, 466, 467, 468, 469, 86, 87, 472, 473, 474, 475, 88, 85, 732, 733, 996, 618, 619, 620, 621, 110, 632, 624, 625, 626, 627, 1012, 622, 623, 1016, 633, 634, 635, 636, 637}, {640, 641, 642, 643, 646, 647, 60, 61, 971, 972, 974, 87, 89, 90, 476, 477, 478, 479, 480, 481, 736, 737, 484, 485, 486, 487, 108, 109, 1015, 628, 629, 1013, 631, 630, 1014, 638, 639}, {640, 641, 1024, 644, 645, 646, 647, 648, 649, 652, 653, 639, 61, 62, 973, 974, 976, 89, 90, 91, 482, 483, 484, 485, 486, 487, 738, 739, 490, 491, 492, 493, 108, 113, 1013, 638, 1023}, {1025, 1026, 644, 645, 648, 649, 650, 651, 652, 653, 654, 655, 767, 660, 661, 62, 63, 975, 976, 978, 90, 91, 92, 488, 489, 490, 491, 492, 493, 496, 497, 498, 499, 113, 114, 766, 1023}, {1025, 1027, 650, 651, 654, 655, 656, 657, 658, 659, 660, 661, 662, 663, 1042, 672, 673, 63, 64, 977, 978, 980, 91, 92, 93, 95, 494, 495, 496, 497, 498, 499, 114, 502, 503, 504, 505}, {64, 65, 512, 66, 513, 514, 515, 1027, 656, 657, 658, 979, 980, 981, 662, 663, 984, 659, 92, 93, 95, 500, 501, 502, 503, 504, 505, 506, 507, 508, 509}, {516, 517, 518, 519, 520, 521, 1031, 774, 524, 525, 526, 527, 775, 1044, 1045, 664, 665, 666, 667, 668, 669, 670, 671, 672, 673, 678, 679, 67, 68, 982, 985, 987, 94, 95, 96, 114, 117}, {512, 513, 514, 515, 1027, 518, 519, 520, 521, 656, 657, 658, 659, 660, 661, 662, 663, 664, 665, 1042, 1044, 668, 669, 672, 673, 66, 67, 983, 984, 985, 92, 93, 94, 95, 114, 510, 511}, {1028, 1031, 1032, 522, 523, 524, 525, 526, 527, 778, 779, 530, 531, 532, 533, 666, 667, 670, 671, 674, 675, 676, 677, 678, 679, 682, 683, 68, 69, 986, 987, 989, 94, 96, 97, 115, 117}, {1028, 1029, 1030, 776, 777, 528, 529, 530, 531, 532, 533, 536, 537, 538, 539, 674, 675, 676, 677, 680, 681, 682, 683, 684, 685, 690, 691, 69, 70, 988, 989, 991, 96, 97, 98, 115, 116}, {1029, 1033, 1043, 534, 535, 536, 537, 538, 539, 542, 543, 544, 545, 680, 681, 684, 685, 686, 687, 688, 689, 690, 691, 692, 693, 70, 71, 990, 991, 97, 98, 99, 993, 112, 116, 762, 763}, {71, 72, 1033, 74, 540, 541, 542, 543, 544, 545, 992, 99, 548, 549, 98, 993, 546, 547, 998, 686, 687, 558, 559, 112, 560, 561, 692, 693, 689, 688, 1021}, {384, 385, 1035, 1039, 552, 553, 42, 43, 554, 555, 564, 565, 694, 695, 568, 569, 696, 697, 73, 75, 994, 995, 100, 997, 111, 379, 756, 757, 374, 375, 376, 377, 378, 1019, 118, 382, 383}, {130, 1037, 1046, 1063, 812, 813, 1073, 698, 570, 571, 700, 701, 702, 703, 704, 705, 578, 579, 580, 581, 699, 572, 573, 574, 714, 715, 76, 77, 575, 712, 713, 1084, 730, 731, 706, 101, 102, 103, 1000, 707, 1002, 106, 122}, {129, 130, 1036, 792, 793, 1070, 1073, 562, 563, 1074, 566, 567, 568, 569, 698, 699, 572, 573, 574, 575, 706, 707, 708, 709, 710, 711, 712, 713, 75, 76, 844, 845, 101, 102, 1001, 1002, 118}, {1037, 1038, 1040, 1046, 700, 701, 702, 703, 576, 577, 578, 579, 580, 581, 584, 585, 586, 715, 587, 714, 77, 78, 590, 81, 591, 80, 594, 595, 598, 599, 596, 597, 101, 103, 106, 1004, 1018}, {1050, 1052, 1053, 806, 807, 730, 716, 717, 718, 719, 720, 721, 82, 83, 722, 723, 600, 601, 602, 603, 604, 605, 729, 728, 608, 609, 610, 611, 731, 104, 105, 106, 1005, 1007, 1009, 121, 122}, {1049, 1050, 1051, 800, 801, 716, 717, 720, 721, 83, 84, 724, 725, 726, 727, 729, 728, 606, 607, 608, 609, 610, 611, 734, 735, 614, 615, 104, 105, 616, 617, 107, 1006, 1007, 1011, 120, 121}, {1040, 1046, 1052, 1063, 700, 701, 702, 703, 704, 705, 714, 715, 718, 719, 80, 592, 82, 593, 596, 597, 594, 595, 722, 723, 602, 603, 604, 605, 730, 731, 101, 103, 104, 106, 1008, 1009, 122}, {1034, 1041, 1048, 1049, 84, 85, 724, 725, 88, 726, 727, 732, 733, 734, 735, 612, 613, 614, 615, 616, 617, 105, 107, 620, 621, 622, 623, 110, 750, 1010, 1011, 751, 120, 634, 635, 636, 637}, {640, 641, 642, 643, 1024, 128, 646, 647, 648, 649, 772, 773, 1054, 1067, 1068, 818, 819, 89, 90, 736, 737, 738, 739, 740, 741, 742, 743, 746, 747, 108, 109, 113, 1013, 1015, 123, 638, 639}, {640, 641, 642, 643, 1054, 1057, 1059, 816, 817, 1017, 87, 89, 736, 737, 740, 741, 744, 745, 746, 747, 108, 109, 110, 748, 749, 752, 753, 1015, 628, 629, 630, 631, 632, 1014, 633, 123, 125}, {637, 1041, 1048, 632, 1057, 802, 803, 1058, 87, 88, 732, 733, 120, 734, 735, 744, 745, 107, 748, 109, 110, 750, 751, 753, 626, 627, 749, 752, 630, 631, 1016, 633, 634, 635, 636, 1017, 125}, {1039, 790, 791, 1047, 1022, 794, 795, 1060, 1062, 552, 553, 554, 555, 556, 557, 694, 695, 696, 697, 73, 100, 759, 111, 112, 754, 755, 756, 757, 118, 119, 760, 758, 761, 1019, 764, 765, 126}, {764, 1033, 765, 784, 785, 1043, 1060, 1061, 550, 551, 554, 555, 556, 557, 686, 558, 688, 559, 560, 561, 692, 693, 687, 689, 690, 691, 73, 74, 98, 99, 126, 111, 112, 754, 755, 116, 760, 761, 762, 763, 1020, 1021, 1022}, {1024, 768, 1026, 769, 644, 645, 646, 647, 648, 649, 774, 775, 652, 653, 654, 655, 788, 789, 128, 1023, 1055, 1056, 770, 771, 1067, 772, 1069, 773, 824, 825, 90, 91, 738, 739, 742, 743, 108, 113, 114, 117, 124, 766, 767}, {768, 1025, 1026, 769, 774, 775, 650, 651, 652, 653, 654, 655, 658, 659, 660, 661, 1042, 1044, 664, 665, 1045, 668, 669, 670, 671, 672, 673, 1055, 91, 92, 94, 95, 113, 114, 117, 766, 767}, {1028, 1030, 776, 777, 1032, 778, 779, 780, 782, 783, 781, 786, 787, 788, 789, 674, 675, 676, 677, 678, 679, 1064, 682, 683, 684, 685, 1066, 1072, 822, 823, 96, 97, 115, 116, 117, 124, 127}, {1029, 1030, 776, 777, 782, 783, 784, 785, 786, 1043, 787, 1061, 680, 681, 682, 683, 684, 685, 1064, 1065, 688, 689, 690, 691, 832, 833, 97, 98, 112, 115, 116, 762, 763, 764, 765, 126, 127}, {768, 769, 770, 771, 774, 1031, 1032, 775, 778, 779, 780, 781, 788, 1045, 789, 666, 667, 668, 669, 670, 671, 1055, 1056, 676, 677, 678, 679, 1066, 94, 96, 113, 114, 115, 117, 124, 766, 767}, {129, 1035, 1036, 1039, 790, 791, 1047, 792, 793, 796, 797, 1070, 1071, 564, 565, 566, 695, 696, 697, 567, 568, 569, 694, 708, 709, 710, 711, 75, 100, 102, 759, 111, 756, 757, 118, 119, 758}, {129, 834, 835, 134, 848, 849, 790, 791, 1047, 792, 794, 795, 793, 796, 797, 798, 799, 1080, 1062, 759, 111, 1071, 756, 757, 758, 119, 760, 118, 761, 126, 1081}, {132, 1048, 1049, 1051, 800, 801, 802, 803, 1058, 804, 805, 808, 809, 1076, 1077, 828, 829, 724, 725, 726, 727, 728, 729, 732, 733, 734, 735, 105, 107, 110, 750, 751, 752, 753, 120, 121, 125}, {132, 133, 1050, 1051, 1053, 800, 801, 804, 805, 806, 807, 808, 809, 810, 811, 814, 815, 1076, 1078, 1079, 716, 717, 720, 721, 722, 723, 726, 727, 728, 729, 860, 861, 104, 105, 120, 121, 122}, {130, 133, 1052, 1053, 806, 807, 1063, 810, 811, 812, 813, 814, 815, 1078, 1084, 1085, 702, 703, 704, 705, 706, 707, 718, 719, 720, 721, 722, 723, 850, 851, 730, 731, 101, 104, 106, 121, 122}, {128, 1089, 1091, 135, 842, 843, 1054, 736, 737, 1059, 740, 741, 742, 743, 746, 747, 108, 109, 748, 749, 816, 817, 1068, 818, 819, 820, 821, 123, 125, 830, 831}, {768, 769, 770, 771, 128, 772, 773, 131, 778, 779, 780, 781, 782, 783, 788, 789, 1056, 1066, 1069, 1072, 1075, 822, 823, 824, 825, 826, 827, 1086, 836, 837, 840, 841, 113, 115, 117, 124, 127}, {132, 135, 1057, 802, 803, 1058, 1059, 805, 804, 816, 817, 820, 1077, 821, 828, 829, 830, 831, 1089, 1090, 862, 863, 744, 745, 746, 747, 748, 109, 110, 749, 752, 753, 751, 750, 120, 123, 125}, {134, 784, 785, 786, 787, 794, 795, 798, 799, 1060, 1061, 1062, 1065, 1081, 832, 833, 834, 835, 1088, 838, 839, 759, 111, 112, 754, 755, 116, 758, 119, 760, 761, 762, 763, 764, 765, 126, 127}, {131, 134, 776, 777, 780, 781, 782, 783, 784, 785, 786, 787, 1064, 1065, 1072, 1075, 822, 823, 826, 827, 1087, 832, 833, 834, 835, 836, 837, 838, 839, 1088, 856, 857, 115, 116, 124, 126, 127}, {128, 770, 771, 772, 773, 131, 135, 1067, 1068, 1069, 818, 819, 820, 821, 824, 825, 826, 827, 1086, 1091, 1094, 840, 841, 842, 843, 858, 859, 738, 739, 740, 741, 742, 743, 108, 113, 123, 124}, {129, 130, 133, 134, 790, 791, 792, 793, 796, 797, 798, 799, 1070, 1071, 1074, 1080, 1082, 1083, 708, 709, 710, 711, 712, 713, 844, 845, 846, 847, 848, 849, 850, 851, 864, 865, 102, 118, 119}, {704, 129, 130, 706, 707, 133, 710, 711, 712, 713, 1082, 844, 845, 846, 847, 850, 851, 705, 101, 102, 122, 812, 813, 814, 815, 1073, 1074, 698, 699, 1084, 1085}, {128, 131, 132, 133, 134, 135, 1075, 822, 823, 824, 825, 826, 827, 1086, 1087, 836, 837, 838, 839, 840, 841, 1092, 1093, 1094, 842, 843, 1095, 852, 853, 854, 855, 856, 857, 858, 859, 860, 861, 862, 863, 864, 865, 124, 127}, {131, 132, 133, 135, 800, 801, 802, 803, 804, 805, 808, 809, 810, 811, 1076, 1077, 1079, 828, 829, 830, 831, 1090, 1092, 1095, 852, 853, 854, 855, 858, 859, 860, 861, 862, 863, 120, 121, 125}, {129, 130, 131, 132, 133, 134, 806, 807, 808, 809, 810, 811, 812, 813, 814, 815, 1078, 1079, 1082, 1083, 1085, 1092, 1093, 844, 845, 846, 847, 848, 849, 850, 851, 852, 853, 854, 855, 856, 857, 860, 861, 864, 865, 121, 122}, {129, 131, 133, 134, 794, 795, 796, 797, 798, 799, 1080, 1081, 1083, 1087, 832, 833, 834, 835, 836, 837, 838, 839, 1088, 1093, 846, 847, 848, 849, 854, 855, 856, 857, 864, 865, 119, 126, 127}, {128, 1089, 1090, 1091, 132, 131, 1094, 135, 840, 841, 842, 843, 1095, 852, 853, 858, 859, 862, 863, 816, 817, 818, 819, 820, 821, 125, 123, 828, 829, 830, 831}]\n"
     ]
    }
   ],
   "source": [
    "blocks = []\n",
    "freedofs = fes.FreeDofs()\n",
    "for v in mesh.vertices:\n",
    "    vdofs = set()\n",
    "    for el in mesh[v].elements:\n",
    "        vdofs |= set(d for d in fes.GetDofNrs(el) if freedofs[d])\n",
    "    blocks.append (vdofs)\n",
    "print (blocks)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`CreateBlockSmoother` can now take these blocks and construct a block Jacobi preconditioner."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "34.84033814053211"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "blockjac = a.mat.CreateBlockSmoother(blocks)\n",
    "\n",
    "lams = EigenValues_Preconditioner(mat=a.mat, pre=blockjac)\n",
    "max(lams)/min(lams)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Multiplicative smoothers and its symmetrized version often yield better condition numbers in practice. We can apply the same code we wrote above for symmetrization (`SymmetricGS`) to the block smoother:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.982606144898276"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "blockgs = SymmetricGS(blockjac)\n",
    "\n",
    "lams = EigenValues_Preconditioner(mat=a.mat, pre=blockgs)\n",
    "max(lams)/min(lams)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Add a coarse grid correction\n",
    "\n",
    "Dependence of the condition number on degrees of freedom can often be reduced by preconditioners that appropriately use a coarse grid correction.  It is also possible to experiment with coarse grid corrections using NGSolve's python interface. We now show how to precondition with a coarse grid correction made using the lowest order subspace of `fes`.\n",
    "\n",
    "In the example below, note that we use `fes.GetDofNrs` again. Previously we used it with argument `el` of type `ElementId`, while now we use it with an argument `v` of type `MeshNode`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0: 00100000000001111111111111111110000000001111111111\n",
      "50: 11111111111111111111111111111111111111111111111111\n",
      "100: 11111111111111111111111111111111111100000000000000\n",
      "150: 00000000000000000000000000000000000000000000000000\n",
      "200: 00000000000000000000000000000000000000000000000000\n",
      "250: 00000000000000000000000000000000000000000000000000\n",
      "300: 00000000000000000000000000000000000000000000000000\n",
      "350: 00000000000000000000000000000000000000000000000000\n",
      "400: 00000000000000000000000000000000000000000000000000\n",
      "450: 00000000000000000000000000000000000000000000000000\n",
      "500: 00000000000000000000000000000000000000000000000000\n",
      "550: 00000000000000000000000000000000000000000000000000\n",
      "600: 00000000000000000000000000000000000000000000000000\n",
      "650: 00000000000000000000000000000000000000000000000000\n",
      "700: 00000000000000000000000000000000000000000000000000\n",
      "750: 00000000000000000000000000000000000000000000000000\n",
      "800: 00000000000000000000000000000000000000000000000000\n",
      "850: 00000000000000000000000000000000000000000000000000\n",
      "900: 00000000000000000000000000000000000000000000000000\n",
      "950: 00000000000000000000000000000000000000000000000000\n",
      "1000: 00000000000000000000000000000000000000000000000000\n",
      "1050: 0000000000000000000000000000000000000000000000\n"
     ]
    }
   ],
   "source": [
    "vertexdofs = BitArray(fes.ndof)\n",
    "vertexdofs[:] = False\n",
    "\n",
    "for v in mesh.vertices:\n",
    "    for d in fes.GetDofNrs(v):\n",
    "        vertexdofs[d] = True\n",
    "        \n",
    "vertexdofs &= fes.FreeDofs()\n",
    "\n",
    "print(vertexdofs)    # bit array, printed 50 chars/line"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Thus we have made a mask `vertexdofs` which reveals all free dofs associated to vertices. If these are labeled $c$ (and the remainder is labeled $f$), then the matrix $A$ can partitioned into \n",
    "\n",
    "$$\n",
    "A = \\left( \\begin{array}{cc} A_{cc} & A_{cf} \\\\ A_{fc} & A_{ff} \\end{array} \\right). \n",
    "$$\n",
    "\n",
    "The matrix `coarsepre` below represents\n",
    "\n",
    "$$\n",
    "\\left( \\begin{array}{cc} A_{cc}^{-1} & 0 \\\\ 0 & 0 \\end{array} \\right). \n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "coarsepre = a.mat.Inverse(vertexdofs)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This matrix can be used for coarse grid correction. \n",
    "\n",
    "*Pitfall!*  Note that `coarsepre` is not appropriate as a preconditioner by itself as it has a large null space. You might get the wrong idea from the results of a Lanczos eigenvalue estimation:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "       1"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "EigenValues_Preconditioner(mat=a.mat, pre=coarsepre)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But this result only gives the Laczos eigenvalue estimates on the *range* of the preconditioner. The preconditioned operator in this case is simply \n",
    "\n",
    "$$\n",
    "\\left( \\begin{array}{cc} A_{cc}^{-1} & 0 \\\\ 0 & 0 \\end{array} \\right)\n",
    "\\left( \\begin{array}{cc} A_{cc} & A_{cf} \\\\ A_{fc} & A_{ff} \\end{array} \\right)\n",
    " = \n",
    " \\left( \\begin{array}{cc} I  & A_{cc}^{-1} A_{cf} \\\\ 0 & 0  \\end{array} \\right),\n",
    "$$\n",
    "\n",
    "which is a projection into the $c$-dofs. Hence its no surprise that Lanczos estimated the eigenvalues of this operator (on its range) to be just one. But this does not imply that the condition number of this preconditioned system is nice.\n",
    "\n",
    "One well-known and correct way to combine the coarse grid correction with one of the previous smoothers is to combine them additively,  to get an **additive two-grid preconditioner** as follows."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "twogrid = coarsepre + blockgs "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This addition of two operators (of type `BaseMatrix`) results in another operator, which is stored as an expression, to be evaluated only when needed.  The 2-grid preconditioner has a very good condition number."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       " 0.993081\n",
       " 0.997768\n",
       " 0.999961\n",
       " 1.34206\n",
       " 1.80314\n",
       " 1.83857\n",
       " 1.96023\n",
       "  1.9828\n",
       " 1.99987"
      ]
     },
     "execution_count": 20,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "EigenValues_Preconditioner(mat=a.mat, pre=twogrid)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
